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We describe a new map-making code for cosmic microwave background (CMB) observations. It 
implements fast algorithms for convolution and transpose convolution of two functions on the sphere 
(Wandelt & Gorski 2001) yj. Our code can account for arbitrary beam asymmetries and can be 
applied to any scanning strategy. We demonstrate the method using simulated time-ordered data for 
three beam models and two scanning patterns, including a coarsened version of the WMAP strategy. 
We quantitatively compare our results with a standard map-making method and demonstrate that 
the true sky is recovered with high accuracy using deconvolution map-making. 



I. INTRODUCTION 

Real microwave telescopes collect distorted informa- 
tion about the cosmic microwave background (CMB) 
anisotropies due to asymmetries in the beam shape |2j 
and stray light from sources such as the Galaxy [J, Q| . 
To correct for these systematic errors we must be able 
to remove the detector response at all orientations of the 
telescope over the whole sky. In an optimal treatment, 
this correction must be applied during the map-making 
step of the CMB data analysis pipeline, before the an- 
gular power spectrum can be reconstructed. The prob- 
lem becomes increasingly important as new generations 
of CMB observations probe for ever fainter signals in the 
CMB sky, and especially as we are preparing to mea- 
sure the polarization of the CMB with high sensitivity. 
We present a complete map-making algorithm, in which 
time-ordered data (TOD) is used to construct a temper- 
ature map and beam distortions are removed. 

We call our approach deconvolution map-making, a 
generalization of existing CMB map-making techniques 
to solve the maximum likelihood map-making problem 
for arbitrary beam shapes. For sufficiently high signal- 
to-noise this technique allows super-resolution imaging 
of the CMB from time-ordered scans. We implement our 
method using the exact algorithms for the convolution 
and transpose convolution of two arbitrary function on 
the sphere - in this case the sky and the beam - as de- 
tailed by Wandelt and Gorski in []J . These fast methods 
for convolution and transpose convolution are efficient 
because they make use of the Fast Fourier Transform 
algorithm. They are guaranteed to work to numerical 
precision for band-limited functions on the sphere. 

Early work on the map-making problem have relied 
on the brute force method of direct matrix inversion. 
However, current and future CMB experiments, like the 
Wilkinson Microwave Anisotropy Probe ( WMAP) |5j 
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and Planck satellite return enormous data sets that 
render the brute force method useless. More recent ad- 
vancements include map-making methods applicable to 
the latest experiments; however, many treat the beam 
like a perfect delta- function {e.g. @,||) or assume a sym- 
metric beam profile {e.g. and thereby relegate the 
problem of treating a non-Gaussian radial response of the 
beam to subsequent stages in the data analysis [lOj. In 
this class, special techniques exist to deal with differen- 
tial measurements like that of the Differential Microwave 
Radiometer (DMR) on the Cosmic Background Explorer 
{COBE) satellite [ill or WMAP A Fourier method 
has been developed 0] to perform deconvolution but 
only for non-rotating asymmetric beams. Lastly, [l4j 
present a method to remove the main beam distortion 
over patches of the sky for asymmetric, rotating beams 
but operate in pixel-space which is computationally more 
expensive than spherical-harmonic-space algorithms for 
the same level of accuracy yj, [l5| . 

We test our algorithm on a simulated foreground- and 
Galaxy-free sky using a standard ACDM power spectrum 
and simulated spherical harmonic multipoles ag m up to 
£ = 128. We also use the first-year WMAP Ka-band tem- 
perature map as our true sky containing Galactic emis- 
sion. 

In section^we present the deconvolution method and 
briefly review a standard map-making method. In sec- 
tion ^Q] we detail the various test cases. Our results for 
the deconvolution method are given, discussed, and com- 
pared with the standard estimates in section I1VI We 
conclude in section Ivl and remark on future directions. 

II. DECONVOLUTION MAP-MAKING 

In order to define our notation we will briefly review 
the path from observations to maps. A microwave tele- 
scope scans the CMB sky according to some scanning 
strategy, effectively convolving the true sky with a beam 
function, and returns a vector, d, containing the Utod 
samples of the time-ordered data. We represent this by 

As = d, (1) 
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where A is the observation matrix, defined below, and s 
is a rip^-vector containing the true sky. 

The matrix A encodes both the scanning strategy and 
the optics of the CMB instrument. Each sample of the 
TOD is modeled as the scalar product of a row of the 
matrix A with the sky s. Each of the utod rows of A 
contains a rotated map of the beam. In a given row the 
beam rotation corresponds to the orientation of the an- 
tenna at the point in time when the sample is taken. We 
will assume the beam shape and pointing of the satellite 
to be known. 

The observation matrix A generalizes the notion of the 
pointing matrix which is often used in expositions of map- 
making algorithms by including both optics and scanning 
strategy. This generalization is necessary for any map- 
making method that accounts for beam functions with 
azimuthal structure. 

The least-squares estimate of the true sky, s, is given 

by 

A T As = A T d. (2) 

The coefficient matrix in this system of equations, A T A, 
is a smoothing matrix and hence ill-conditioned. Invert- 
ing it to solve Eq. J2J therefore poses a problem. 

We describe here a regularization technique for dealing 
with this problem. We split off the ill-conditioned part 
of A by factoring the convolution operator into A = 
BG where G is a simple Gaussian smoothing matrix, 
represented in harmonic-space by 

G « = e xp (^f±I>), (3) 

where a = FWHM/V81n2. 

Substituting the factorization into Eq. J5J), we get 

G T B T BGs = G T B T d (4) 
B T Bx = B T d (5) 

where we are solving for x = Gs so as not to reconstruct 
the sky at higher resolution than that of the instrument. 

Equation J3J) is exact if the noise is stationary and un- 
corrected in the time-ordered domain. For a more gen- 
eral noise covariance matrix in the time-ordered domain, 
N, the normal equation is modified as follows 

A T N _1 As = A T N _1 d. (6) 

We proceed, considering only white noise in this paper; 
however, it is straightforward to generalize to non-white 
noise as indicated in Eq. J§}. Indeed, the matrix- vector 
operations required for this generalization have already 
been implemented in publically available map-making 
codes (e.g. MADMAP [l6|). 

A. Fast Convolution on the Sphere 

We now briefly review the relevant formalism for gen- 
eral convolutions on the sphere and refer the reader to 



for the full details of the fast convolution and transpose 
convolution of two functions in the spherical harmonic 
basis. The convolution of a band-limited beam function 
6(7) with the sky s(t*) is given by the following integral 
over all solid angles 

T(^ 2 ,Q,^ 1 ) = J dShtlDfa.QMMWstf) (7) 

where D is the operator of finite rotations 23] . In spher- 
ical harmonic space this becomes 

Tmm'm" = ^ si m d e mm (8 E )d l m , m „{0)b* im n . (8) 

i 

Analogously, the transpose convolution of T($ 2 , 6, $1) 
is given by 

y*(7)= / d^dQd^bi^Q^^b^YT^Q,^), 

(9) 

and in spherical harmonics 

ytm = X! dmrn' (QE)d e m , m „ (8)be m „T mm > m » . (10) 
m'm" 

An important feature of our approach is that it econ- 
omizes the computational effort if the beam is nearly az- 
imuthally symmetric. The parameter of the method that 
sets the degree to which asymmetries of the beam are 
taken into account is TO max , the maximum m" in equa- 
tions (fTT7f> and © ■ For m max = we recover the compu- 
tational cost of simple spherical harmonics transforms, 
C(^max)- Since TO max is bounded from above by £ max , 
the computational cost of the method never scales worse 
than C(^ ax ). For a mildy elliptical beam, we antici- 
pate that just including the m max = and m max = 2 
terms will suffice, since the m max = 1 term vanishes by 
symmetry. 

For clarity, we now rewrite Eq. (JSJ) in the compact 
spherical-harmonic basis (summing over repeated in- 
dices) 

-^■L'At'mm'm"^mm'm"LM x LM — ^L' M' mm' m'^mm' m" j 

(11) 

where A T acting on T mm ' m " is given by Eq. (|10|) and B 
acting on xlm is given by Eq. (JSJl. 

To make matters even more concrete, we now explicitly 
describe the steps required to simulate time-ordered data 
d from a map ("simulation"). We convolve the beam bg m 
with the map ai m to obtain T mm t m ». Then we inverse 
Fourier transform the T mm / m » to get T(<i>2, @, ^i)- Next, 
we must account for the scan path (<I>2(i), O(t), $i(i)), 
where $2 and specify the position on the sphere and 
$1 specifies the orientation of the beam. This is achieved 
by extracting those values in T($2, O, ^l) which fall on 
the scan path whenever the instrument samples the sky. 

As a second example we describe how to compute the 
right hand side of Eq. JSJl. Start with the TOD d. For 
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each sample in d, the scanning strategy specifies the ori- 
entation ($2, ©) 3>i)< The sampled temperature is added 
into the element of an initially empty array which is 
identical in size and shape to the array which stored 
T($ 2 ,e,$i). We have effectively binned the TOD d, 
according the position and orientation of the beam on 
the sky. Let us therefore refer to this operation as "bin- 
ning" . In order to minimize discreteness effects due to 
the gridded representation of T(<E> 2 , 0, $i), more sophis- 
ticated interpolation techniques could be implemented. 
Additionally, the resolution of the grid into which the 
data is binned may be increased. 



B. Solving the Deconvolution Equations 

To obtain the optimal map estimate we numerically 
solve the linear system of equations in Eq. (0 for xg m . 
We have a choice between direct and iterative solution 
methods. An iterative method is advantageous compared 
to a direct method (such as Cholesky inversion) if the cost 
per iteration times the number of iterations required to 
converge to sufficient accuracy is less than the cost of the 
direct method. 

For the problem sizes of current and upcoming CMB 
missions, where the map contains a number of pixels 
Upix ~ 10 6 -10 7 direct solution methods would be pro- 
hibitive for two reasons. Firstly, the required number of 
floating point operations scales as v? vix . Secondly, the 
amount of space required to store the coefficient matrix 
and its inverse scales as nz ix . Therefore direct solution 
exceeds the capabilities of modern supercomputers by 
several orders of magnitude. For the Planck mission di- 
rect solution would require of order 10 21 floating point 
operations and hundreds of Terabytes of random access 
memory. 

We therefore advocate using an iterative technique, the 
Conjugate Gradient (CG) method [HI- The CG method 
is well suited to this problem. It solves linear systems 
with symmetric positive definite coefficient matrix and 
has advantageous convergence properties compared to 
other iterative methods such as Jacobi method |l7j. In 
order to apply the CG method we must be able to apply 
the coefficient matrix on the left hand side of Eq. (J2J) to 
our current guess of the solution x. In order to do so we 
simply perform the two operations of "simulation" and 
"binning" in succession. The fast convolution and trans- 
pose convolution algorithms allow computing the action 
of the coefficient matrix on a map without ever having 
to store the matrix coefficient in memory. 

It is desirable to minimize the number of iterations the 
CG method requires to converge to a given level of accu- 
racy. This can be done by "preconditioning" the system 
of equations. Preconditioning amounts to multiplying on 
both sides with an approximation of the inverse of the 
coefficient matrix and solving this modified system. As 
long as the preconditioner is non-singular the solution 
will be the same for the original and the preconditioned 



systems, but for a well-chosen preconditioner the num- 
ber of iterations can be reduced significantly. A natural 
choice of the preconditioning matrix which we used to 
obtain the results in this paper, is the diagonal matrix 
[diag(A T A)]~ 1 , which, for a delta-function beam, is just 
the inverse of the number of hits per pixel. 

At every iteration we have an approximate solution x 
of Eq. We assess convergence by computing the ratio 
of L 2 norms 

£ 2 [B T Bx B T d] 

£ 2 [B T d] 1 ' 

where £ 2 [x] = yjx • x|. 



C. Standard Map-Making: Brief Review and 
Critique 

In order to compare our results to traditional tech- 
niques we also implemented a traditional map-making 
code that solves the normal equation (Eq. iffll) assuming 
an azimuthally symmetric beam. In this implementation 
the observation matrix A becomes the pointing matrix, 
containing only a single entry on each row corresponding 
to the direction in which the main beam lobe is pointing 
at the time of sampling. Standard map-making therefore 
reconstructs a map which is smoothed by an effective 
beam whose shape varies as function of position on the 
map. This variation depends on the scanning strategy. 
More precisely, at any given position on the estimated 
map the effective beam shape depends on the various ori- 
entations of the beam as it passed through this position 
during the scan. 

For uncorrelated noise and an azimuthally symmetric 
beam the solution of the normal equation is simple to 
compute: bin the TOD into discrete sky pixels, sum- 
ming over repeated hits, and dividing through by the 
number of hits per pixel. Numerical implementations of 
this algorithm and its generalization to correlated noise 
have been described in the literature However, 
all of these treatments assume azimuthally symmetric 
beams. For experiments with highly asymmetric beams 
and where contamination from the Galaxy is picked up 
in the sidelobes, we expect that this method will not 
fare well against our deconvolution method which also 
removes artifacts due to these optical systematics. We 
use the same TOD and scan path as for our deconvo- 
lution method. Here, the data is binned into pixelized 
maps, rather than into the T(<i>2, O, $1) grid. Unless oth- 
erwise stated we use the HEALPix pixelization scheme 
|18| with resolution parameter nside = 64. The angular 
scale of a pixel is therefore just under 1°. Recall that our 
rcgularization method returns a smoothed map with an 
effective, azimuthally symmetric Gaussian beam. Thus, 
in order to compare the two methods we must make a 
similar modification to our standard map-making. We 
read out the resulting ag m of our standard map (using 
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the HEALPix anafast routine), after the binning step, 
and modify them in the following way 



where Bg is the beam power spectrum and Gi is given in 
Eq. ©. 

III. TEST CASES 

In this section we detail our tests and comparisons 
of the deconvolution and standard map-making meth- 
ods. For the purposes of testing our method, we create 
several mock beam models bg m . We test three possible 
beam shapes which break azimuthal symmetry progres- 
sively strongly, two scanning patterns, and skies with and 
without Galactic emission. 

The first beam is a simple model of a sidelobe; it is 
composed of a Gaussian beam of FWHM = 1800' rotated 
at 90° to another Gaussian beam of FWHM = 180'. Both 
the main beam and the sidelobe are normalized such that 
they integrate to one. The second beam models a (some- 
what exaggerated) elliptical shape, composed of two iden- 
tical Gaussian beams with FWHM = 180' whose centers 
are on both sides of the optical axis, separated by 180'. 
The third beam is composed of two identical Gaussian 
beams (FWHM = 180') rotated at 140° from each other; 
we refer to this as the two-beam model. This case is 
motivated by the design of the WMAP observatory . 

We set the asymmetry parameter m max for our three 
cases (sidelobe, elliptical, and two-beam) to 8, 38, and 
128, respectively. 

Following we first considered a basic scan path 
(BSP) in which the beam scans the full sky on rings of 
constant longitude with no rotation about its outward 
axis. To be clear, for the case of the sidelobe beam, the 
smaller beam follows this ringed-scan while the larger 
beam remains fixed at the equatorial longitude. Simi- 
larly, in the two-beam model, one beam follows the ring- 
scan while the other rotates in smaller circles 140° away. 
The central lobe therefore covers the whole sky, while the 
offset beam remains within a band of ±50° centered on 
the ecliptic plane. The elliptical beam simply follows the 
ring-scan, and is oriented such that its long axis remains 
perpendicular to the lines of longitude. 

A more realistic observational strategy has a beam 
that revisits locations on the sky in different orienta- 
tions. Therefore, we model the one-year WMAP scan 
path followed by one horn. The WMAP scan strategy 
also covers the full sky and includes a spin modulation 
of 0.464 revolutions per minute and a spin precession of 
one revolution per hour [l^. We used a scaled-down 
model of the WMAP scan in which the spin modula- 
tion is 0.00232 revolutions per minute and a step size of 
about 46 seconds (roughly 562 samples per period). This 
produces a pattern very similar to the spirograph-type 
pattern shown in Fig. 4 of 0|. We refer to this as the 



WMAP-like scan path (WSP). The WSP has about six 
times as many samples as the BSP. For this strategy, the 
spirograph pattern is followed by the small beam of the 
sidelobe, the elliptical beam, and both beams of the two- 
beam model. In the two-beam case both beams are offset 
from the spin axis of the satellite, to mimic the WMAP 
scanning geometry. It is not differential in nature, since 
both beams have positive weight. 

We test each beam (sidelobe, elliptical, and two-beam) 
with both scanning patterns (BSP and WSP) on a sky 
without Galactic emission. We refer to these as the six 
main test cases. 

In reality, CMB experiments will also pick up signal 
from the Galaxy. We use the first-year WMAP Ka- 
band temperature map, degraded to an nside of 64 and 
smoothed with a Gaussian beam of FWHM = 180' as 
our model of the true sky with Galactic emission. For 
our last test case, we convolve it with the sidelobe beam. 

For each test case, we assume that the beamshapes 
of the instrument are known and use the deconvolution 
method to deconvolve the map with the same beam that 
the true sky was originally convolved with. We attempt 
to recover features in the map corresponding to the small- 
est scale features of our test beams. We therefore set the 
width of our regularization kernel, represented by the 
matrix G in Eqs. © and (@J, to FWHM = 180' in every 
case. We compare our map estimates to the true sky, 
smoothed by the regularization kernel. When we refer 
to the "true" sky in the following we mean this kernel- 
smoothed input sky. 



IV. RESULTS AND DISCUSSION 

We present the results of the deconvolution algorithm 
for the six main test cases in the form of residual maps. 
We compare these results to the results from standard 
map-making by examining their power spectra and by 
calculating the root mean square (RMS) difference be- 
tween the estimated and true sky. For the tests that 
include the Galactic signal we show the actual map esti- 
mates. 

In Fig. we plot ratios of the power spectra of the 
residual maps (both standard and deconvolved) and the 
power spectrum of the input map. The BSP (WSP) re- 
sults are plotted in the left (right) column. The solid 
(dashed) lines represent the relative difference in Cg be- 
tween the deconvolved (standard) map and true sky map. 
The standard map- making algorithm failed to give mean- 
ingful results for the two-beam test. We therefore ex- 
cluded this case from the plot. 

For all cases we chose to present the results after a 
fixed number of iterations to show the impact of scanning 
strategy and beam pattern on the condition number of 
the map-making equations. We find that the deconvo- 
lution algorithm outperforms standard map-making by 
orders of magnitude in accuracy. 

For a fixed number of iterations, the BSP tests per- 
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formed less well than the WSP tests. The two-beam BSP 
and, to a lesser extent, the elliptical beam BSP test cases 
have not converged to sufficient accuracy. 

There are several possible causes for this behaviour. 
The BSP leads to an extremely non-uniform sky cover- 
age. Also, the BSP visits each pixel in a narrow range 
of beam orientations. Further, the number of sky sam- 
ples is smaller for the BSP case than for the WSP case 
(as noted in section ITTTf) . All of these aspects can con- 
tribute to increasing the condition number of the normal 
equation, which in turn leads to smaller error decay per 
iteration of the preconditioned CG solver. 

In Table 1 we summarize the RMS difference between 
the reconstructed and true sky. The RMS values are 
computed using the standard deviations the residual and 
true maps: 



RMS 



stdev(residual map) 
stdev(true map) 



where 



residual map = estimated map — true map. 



(14) 



(15) 



The RMS values reflect the trends seen in the spectra 
in Fig. n The residual maps are shown in Fig. In 
order that the scale of the axes on the residual maps are 
meaningful, we also show the true sky map. 



Beam 


BSP 


WSP 


Standard 


Deconvolved 


Standard 


Deconvolved 


Sidelobe 

Elliptical 

Two-beam 


0.257467 
0.186262 
N/A 


0.000111903 

0.0207830 

0.102778 


0.178828 
0.129715 

N/A 


3.13741e-07 
2.25602e-05 
1.08579e-06 



TABLE I: Fractional RMS error for each of the six main test 
cases. 

Achieving a stably converging iterative solution 
method for the deconvolution problem is a success of our 
rcgularization technique. The convergence of our itera- 
tive solver as function of iteration number is plotted in 
Fig. El 

In order to be able to compare the performance of our 
method for different beam patterns and scanning strate- 
gies we make the deliberate choice of limiting the number 
of iterations to 100 and comparing the best results ob- 
tained up to this point. Since our error estimate contin- 
ues to drop stably (except in two cases, where we reach 
the single precision numerical accuracy floor after ~ 25 
and ~ 70 iterations) it is clear that the accuracy of the 
reconstruction can be improved by allowing the system 
to iterate further, or by choosing a more sophisticated 
pre-conditioner . 

Our final test case consists of a model sky with Galaxy 
emission convolved with the sidelobe beam over the WSP. 
We present the output maps of both the standard and de- 
convolution methods in Fig. 0] where the maps are shown 
in ecliptic coordinates. One can see that the standard 



map contains a distorted image of the Galaxy and that 
the deconvolved map is virtually identical to the true 
map. 



V. CONCLUSIONS AND FUTURE WORK 

We have presented a deconvolution map-making 
method for data from scanning CMB telescopes. Our 
methods remove artifacts due to beam asymmetries and 
far sidelobes. We compare our technique with the stan- 
dard map-making method and demonstrate that the true 
sky is recovered with greatly enhanced accuracy via the 
deconvolution method. Deconvolution map-making re- 
covers features of the CMB sky on the smallest scale of 
the beam, thereby achieving a form of super-resolution 
imaging. This extracts more of the information content 
in CMB data sets. 

One of the key difficulties encountered in deconvolution 
problems is that the systems of linear equations we need 
to solve are very nearly singular. We solve this problem 
by introducing a regularization method which allows us 
to solve the systems stably and recover maps at a uniform 
resolution and with an effective beam that is azimuthally 
symmetric and has a Gaussian profile. 

We tested the convergence speed of two particular 
scanning strategies and found that the WMAP-like scan 
is superior to the basic scan in both rate of convergence 
and true-sky recovery. We hypothesize that this is due to 
the nature of the BSP, where the poles are the location 
of the only beam crossings and receive many more hits 
than the rest of the sky. In addition, our implementation 
of the BSP had a smaller number of samples overall than 
our implementation of the WSP. 

We have also shown the relevance of this algorithm to 
the WMAP mission by demonstrating its operation using 
a WMAP-like scanning strategy and a two-beam model 
which, while not differential, resembles the telescope ori- 
entations of the WMAP spacecraft. Our results underline 
the qualities of the WMAP scanning strategy compared 
to a BSP strategy for deconvolution map-making. 

In order to decouple from issues that are not directly 
related to the optical performance of CMB instruments 
we did not consider the effects of noise in our simulations. 
For a realistic assessment of the performance of our meth- 
ods on real data this needs to be added. In particular, the 
choice of scale for the regularization kernel will depend 
on weighing the benefits of increased resolution against 
increased high-frequency. 

Recently several groups have published CMB polar- 
ization results on the EE power spectrum [20|. WMAP 
released an all-sky analysis of the TE cross-correlation in 
the first year data j^. We eagerly anticipate the large- 
angle polarization data from WMAP in the impending 
release of the second year data. In a few years' time, po- 
larimeters on board the Planck satellite will collect data. 
Owing to the difficulty of separating the two polariza- 
tion modes, we expect that polarimetry experiments will 
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Sidelobe QSP 



Sidelobe WSP 




FIG. 1: Ratios of the spectra of the residual map to the spectrum of the input map for each of the beam models and both 
scanning strategies. The BSP results are plotted in the left column and the WSP results are plotted in the right column. 
Results for the sidelobe, elliptical and two-beam beam are shown in the top, middle, and bottom panels, respectively. The solid 
lines correspond to deconvolved spectra and the dashed lines correspond to the standard spectra. 



be very sensitive to beam asymmetries and stray light. 
Future measurements of the tiny B-mode polarization 
will require both exquisite instruments and sophisticated 
analysis tools. We hope that the generalization of our 
methods to polarized map-making will be useful for mak- 
ing maps of the polarized microwave sky. It has already 
been shown [2^ that little modification to the Wandelt- 
Gorski method of fast all-sky convolution is needed to 
accomodate polarization data. 
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all cases after a fixed number of iterations to show the impact of scanning strategy and beam pattern on the condition number 
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FIG. 3: Convergence rates of the preconditioned conjugate gradient solver for each test case. The left panel refers to the basic 
scan path and the right panel to the WMAP-like scan path. The solid lines correspond to the sidelobe beam, dotted lines to 
the elliptical beam, and dot-dashed to the two-beam model. 
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FIG. 4: Deconvolving the effects of a large sidelobe in simulated observations of the WMAP Ka band map, using the coarsened 
WMAP scanning strategy described in the text. The top map is the input sky map, the middle map is the standard map-making 
result, and the bottom map is the deconvolved result. 



